function dataout = rmoutliers_sd(data, dimension, sdthresh, nmin)
% rmoutliers_sd(data, dimension, sd, nmin)
%   data -- 2D matrix points x sweeps, or sweeps x points
%   dimension -- dimension from which to remove elements (usually the sweep
%       dimension which is 2 for points x sweeps, or 1 for sweeps x points)
%   sdthresh -- Standard deviation threshold for outliers
%   nmin -- Number of outlier points required for element to be removed.
%
% Removes outlier sweeps from a 2D matrix of data.  Sweeps are removed if
% there are enough points in the sweep that are outside a standard
% deviation threshold.  Standard deviation is calcuated at each time point
% separately.

% Subtract sweep triggered mean from each sweep
normdata = data - mean(data, dimension);

% Threshold based on the standard deviation across sweeps at each time point.
thresh = sdthresh * std(normdata, 1, dimension);

% Outliers are above thresh or below -thresh.  Sum occurances accross the
% opposing dimension to get a sum for each sweep.
sumoutliers = sum((normdata > thresh) | (normdata < -thresh), mod(dimension,2)+1);

% Keep all sweeps that are below the mimimum required count
keepindex = sumoutliers < nmin;

% Return selected sweeps.
if dimension == 1
    dataout = data(keepindex, :);
else
    dataout = data(:, keepindex);
end